Tests for one dimensional advection-diffusion are failing on GPUs. I have used similar code to the two dimensional tests which pass so not sure exactly why. Looks to be the same problem for all four tests, a method issue with mul! in set_c! when using dev = GPU

MethodError: no method matching mul!(::CuArray{ComplexF64, 1, CUDA.Mem.DeviceBuffer}, ::CUDA.CUFFT.rCuFFTPlan{Float64, -1, false, 1}, ::Vector{Float64}, ::Bool, ::Bool)

@navidcy I am still quite unfamiliar with the GPU programming in julia and have not been able to fix the 1D tests I added on the GPU, though they all work and pass on the CPU. I thought initially that somehow the array type of c0 was causing the issue (it is a Vector not a CuArray) but after initialising it is a CuArray I still got an error. If you have seen something like this before could you let me know where the error is and I will do my best to fix it! No rush on this!!

One way I thought Problem could use multiple dispatch for the 1D, 2D and (soon) 3D problem is to make the advecting flow an argument instead of a kwarg. In this way it might look something like

function Problem(dev, advecting_flow::Function; kwargs...)
  # Set up one dimensional advection-diffusion `Problem`
end

function Problem(dev, advecting_flow::Vector{Function}; kwargs...)
  if length(advecting_flow) == 2
   # Set up two dimensional advection-diffusion `Problem`
  elseif  length(advecting_flow) == 3
   # Set up three dimensional advection-diffusion `Problem`
  else
   # throw error
  end
end

So we find the dimension of the problem from the advecting_flow argument (though this means the user must pass in this argument and it must be a function).

There is probably a neater/simpler way but this was what I came up with..

I'll have a look at the PR soon!

I'll have a look at the PR soon!

No rush!

Tests all pass now by using x = ArrayType(dev)([gr.x[i] for i ∈ 1:gr.nx]). This is fine though I think it would be better to add

function gridpoints(grid::OneDGrid{T, A}) where {T, A}
  X = [ grid.x[i] for i=1:grid.nx ]
  
  return A(X)
end

to FourierFlows.jl as then x = gridpoints(grid) would give the desired array for a OneDGrid and match the usage of the higher dimensions. I could open a PR to do this if you think it is worth it?

I have changed the way the methods are called for TracerAdvectionDiffusion.Problem. The Problem function now takes an advecting_flow argument rather than u and v as parameters. This is a breaking change. The 1D method now takes a single function for the advecting_flow,

function Problem(dev, advecting_flow::Function;
    nx = 128,
    Lx = 2π,
     κ = 0.1,
    dt = 0.01,
stepper = "RK4",
steadyflow = false,
     T = Float64
)
# Setup OneDProblem
end

For the two dimensional Problem a NamedTuple with the components of the advecting_flow is passed in

function Problem(dev, advecting_flow::NamedTuple=(u=noflow, v=noflow);
          nx = 128,
          Lx = 2π,
          ny = nx,
          Ly = Lx,
           κ = 0.1,
           η = κ,
          dt = 0.01,
     stepper = "RK4",
  steadyflow = false,
           T = Float64
  )
# Set up TwoDProblem
end

It could be extended to include a Problem in 3 dimensions by

function Problem(dev, advecting_flow::NamedTuple; parameters)

 if length(advecting_flow)== 2
  #Setup two dimensional `Problem`
 elseif length(advecting_flow ) == 3
  #Setup three dimensional `Problem`
 else
  #Throw error
 end
end

Not sure what is best or if you want to implement this breaking change!

Here is one way we could structure things with an abstract type for the AdvectingFlow. In the TracerAdvectionDiffusion module we could add

"Abstract supertype for an advecting flow"
abstract type AdvectingFlow end

The the struct for the OneDAdvectingFlow could look like

struct OneDAdvectingFlow <: AdvectingFlow
  "Function for the x-component of the advecting flow"
           u :: Function 
  "Whether of not the flow is time dependent"
  steadyflow :: Bool
end

with a constructor that is exported and has default values like is currently set in Problem,

OneDAdvectingFlow(; u = noflow, steadyflow = false) = OneDAdvectingFlow(u, steadyflow)

Then when setting up a TracerAdvectionDiffusion.Problem if the defaults wanted to be used we could (I have not coded this up and tried) do

# set up params
# use default values
advecting_flow = OneDAdvectingFlow()
prob = TracerAdvectionDiffusion.Problem(dev, advecting_flow; params)

or choose some other function for the flow

# set up params
# set a time dependent flow
u(x, t) = #some function
advecting_flow = OneDAdvectingFlow(; u = u, steadyflow = false)
prob = TracerAdvectionDiffusion.Problem(dev, advecting_flow; params)

Inside the module this would just mean changing u -> advecting_flow.u and steadyflow -> advecting_flow.steadyflow or creating new variables u = advecting_flow.u and steady = advecting_flow.steadyflow.

This could generalise to the 2D and 3D cases quite easily and might be less of a rewrite of the code than using a grid for the different methods of Problem. I am not set on this way, have just been thinking about it and this is what I came up with!

@navidcy

I have updated the TracerAdvectionDiffusion module to use the AdvectingFlow abstract type for the methods of Problem. It is implemented as outlined in #55 (comment) above. I think it is a reasonable structure for the module and setting up a Problem. It is certainly better than what I had previously!

sorry for slacking on this...

I liked the One/TwoDAdvectingFlow!

sorry for slacking on this...

No worries!

Nice work! Thanks!

feel free to merge when tests pass

and tag a new release, you can do that by going to the merge commit and just commenting @juliaregistrator register

Comments on pull requests will not trigger Registrator, as it is disabled. Please try commenting on a commit or issue.

haha, just me saying that got the tagbot triggered :)

feel free to merge when tests pass

Will do! Will look into adding an option for using a ThreeDAdvectingFlow soon and maybe #2 if I have the time!